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Abstract 

In the Hartree-Fock approximation the PauH exclusion principle leads to a Schrodinger Eq. of 
an integro-differential form. We describe a new spectral noniterative method (S-IEM), previously 
developed for solving the Lippman-Schwinger integral equation with local potentials, which has 
now been extended so as to include the exchange nonlocality. We apply it to the restricted case 
of electron-Hydrogen scattering in which the bound electron remains in the ground state and the 
incident electron has zero angular momentum, and we compare the acuracy and economy of the new 
method to three other methods. One is a non-iterative solution (NIEM) of the integral equation 
as described by Sams and Kouri in 1969. Another is an iterative method introduced by Kim and 
Udagawa in 1990 for nuclear physics applications, which makes an expansion of the solution into 
an especially favorable basis obtained by a method of moments. The third one is based on the 
Singular Value Decomposition of the exchange term followed by iterations over the remainder. The 
S-IEM method turns out to be more accurate by many orders of magnitude than any of the other 
three methods described above for the same number of mesh points. 
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I. INTRODUCTION. 

A new very efficient and stable method for solving the Schrodinger equation has recently 
been developed [Q. This method (S-IEM) solves the Lippman-Schwinger integral equation 
rather than the differential Schrodinger equation because the numerical errors for the solution 
of the former are inherently smaller than for the latter, and by making a spectral expansion of 
the solution into Chebyshev polynomials it acquires additional excellent accuracy properties. 
This method also avoids the usual drawback of numerical solutions of integral equations, 
namely, the need to invert large non-sparse matrices which represent the discretized form 
of the equation. It achieves this by dividing the integration interval into partitions and by 
making use of the semi-separable nature of the Green's function in configuration space. The 
accuracy of the S-IEM has been tested by comparing it to the solution of the differential 
Schrodinger equation for various cases, such as the scattering of cold atoms 0, and for 
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tunneling through a barrier , but not by comparing it to existing solutions of the integral 
Lippman-Schwinger equation. The purpose of the present paper is two-fold: a) to extend the 
S-IEM to the case that there are non-local exchange terms present. This is possible 0], 
without losing the original accuracy and stability features since the kernel of the exchange 
integral is semi-separable (see below). And b), to compare the accuracy of the extended 
S-IEM to conventional methods of solving an integral equation in configuration space. 

The scattering of an electron from a hydrogen atom offers a good opportunity for per- 
forming such a test, because, integral equation methods have been developed in the past in 
order to include the integral exchange terms required by the Pauli exclusion principle, as 
is described further below. The present test calculation is physically not realistic since it 
does not allow for the polarization of the bound electron cloud by the incident electron, or 
the ionization of the atom. But it is still sufficiently close to realistic so as to serve as an 
adequate test for the comparison of different algorithms. Realistic calculations which allow 
for the polarization of the electron cloud involve coupling between many channels 0, or 
else the solution of a two-dimensional differential equation |^, which is much beyond the 
scope of the present work, and is not necessary in order to demonstrate the power of our 
new method. 

The three methods for solving the non-local Schrodinger equation in the presence of the 
exchange terms we compare our S-IEM with are as follows. One of these methods |^ solves 
an integral equation which is very similar to the one we solve, with the main difference that 
it uses a trapezium type method for numerically expressing the integrals, rather than the 
spectral Chebyshev method used in the S-IEM. Another is an iterative method introduced 
by Kim and Udagawa in 1990 for nuclear physics applications, which makes an expansion 
of the solution into an especially favorable basis obtained by a method of moments. The 
third is a variation of a conventional iterative method for including the exchange term. It 
consists in expanding the exchange kernel into a small number of separable terms by means 
of the Singular Value Decomposition and then iterating over the remainder |]10[. Since there 
is a rich literature on methods developed for taking the exchange terms into account, a 
review of some of the methods most relevant to the S-IEM will be described below. 

The scattering of electrons from atoms or molecules has been the subject of investigation 
ever since quantum mechanics was introduced, and the research continues unabated, mainly 
concerning the scattering of polarized electrons [|r^ or of high energy photons from 
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atoms or molecules, or charge transfer in ion-molecule scattering [O] . The latter can also be 



calculated by using advanced three-body formalisms |T^. One of the features which gives 
computational difficulties in solving the corresponding Schrodinger equation is the Pauli 
exclusion principle which requires that the wave function of the incident electron be anti- 
symmetric with the wave functions of the electrons in the target atom or molecule. In the 
Hartreee-Fock formulation this requirement leads to the presence of non local terms in the 
(coupled) differential equations of the form 

K{r, r') u{r') dr', (1) 





where u{r) is the required solution and K is the integration kernel due to exchange. In 
the early investigations |T5[ the exchange terms were taken into account iteratively, by 
using Green's functions defined by the local part of the potential. However, under certain 
conditions these iterations do not converge |]l5l , |]l6l . Methods to accelerate the convergence 



have been introduced, but such methods tend to be cumbersome and unpredictable. An 
improved iteration procedure can be obtained by means of a separable representation of the 
integral kernel 0], as is described further below, but this method has its limitations as well. 
An interesting method to include the exchange terms non-iteratively by approximating them 



in terms of a separable representations has been developed by Schneider and Collins [17 



Methods to solve the non-local Schrodinger equation non iteratively and rigorously have 



also been developed [18|. One of the oldest ones originated with I. Percival and R. Marriott 
[T9| . It consists in constructing auxiliary functions which are solutions of the differential 
equation in the presence of several types of inhomogeneous terms and then constructing 
the exact solution (which take into account the integral exchange terms) by means of linear 
combinations of the auxiliary functions. A variation of this method was developed by Lamkin 



and Temkin [ffBl for their Polarized Orbitals procedure. In this method the semi-separable 



nature of the exchange kernel K 

K{r,r') = A{r) B{r') for r' < r 

K{r,r') = A{r')B{r) for r' > r (2) 
is exploited. The integral of Eq. (|1]) then becomes 

oo pr poo 

K{r,r')u{r')dr' = A{r) / B{r') u{r') dr' + B{r) / A{r') u{r') dr' (3) 

Jo Jr 
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which can be rewritten as 

K{r,r')u{r')dr' = A{r) I B{r') u{r') dr' - B{r) I A{r')u{r') dr' + CB{r). (4) 



Again auxihary functions can be obtained by solving the differential integral equation with 
the constant C set equal to zero or unity. The true solution u{r) is obtained by means of a 
linear combination of the auxiliary functions, and the constant C = A{r') u{r') dr' as well 
as the coefficients of the linear combination can be obtained via the solution of an algebraic 
equation. The advantage of this method ]16|] is that, when the integral over the kernel goes 
from to r, as is the case in Eq. @, the solution of the integro-differential equation can 
be performed easily by starting at the origin and increasing the upper limit gradually from 
one meshpoint to the next. With a step size of 0.05 the authors obtain accuracies to about 
three significant figures with this method. 

A method which yields five significant figures with the same step size of 0.05 is obtained by 



Kouri and co-workers @|. The main difference from Temkin's method [|T6| is that the authors 
first transform the integro-differential equation into a Lippman-Schwinger integral equation, 
because integral equations have greater numerical stability than differential equations. They 
again transform the integrals, which originally extend from to oo into integrals from 
to r, plus inhomogeneous terms, thus obtaining Volterra integral equations of the second 
kind. They then easily obtain auxiliary solutions to auxiliary Volterra equations by stepping 
progressively from the origin to increasing values of r, similarly to what is done in the method 
of Temkin .The exact solution, and the respective constants, can then be determined in a 



algebraic way similar to what is done in Ref [|T6l . Smith and Henry have also developed 
methods to solve the Volterra type integral equations non iteratively. These methods are 
generally called NIEM, where the "N" stands for "Non-iterative". Collins and Schneider 
also have examined the NIEM form of the non-local Schrodinger equation, but without 
transforming it into a Volterra type. By discretizing the integral via the trapezoidal rule, 
they obtain a linear algebraic equation for the wave functions at the mesh points, and for 
this reason the method is called (LA). The formulation of the initial equation to be solved 
by our new method is very similar to that of [0. The main differences, to be discussed 
below, arise from the numerical techniques used in the solution of these equations (NIEM). 
For the case of the exchange terms which are due to the Coulomb interaction, as is the case 
for most atomic physics calsulations, the exchange terms can also be replaced by coupling 
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to a set of "pseudo-states, a technique which is made use of in the work of Weatherford, 
Onda and Temkin |^T|, and is also illustrated in the present work. 

Other methods have been developed to solve the electron-atom scattering equations. One 
consist in introducing a set of basis functions such as Laguerre polynomials, and expanding 
both the solution and the target states into this basis [E^. Another such expansion basis 
utilizes sturmian functions In these procedures the exchange integrals can be carried 
out, since they contain the known basis wave functions. Finite element representations have 



also been applied |2J]. The R-matrix approach is also well developed . 

As already briefly mentioned above, our new "spectral integral equation method", S- 
lEM, transforms the differential equation into an equivalent Lippmann-Schwinger integral 
equation through the use of Green's functions, similar to what is done in the older approaches 
described above. It differs from the older methods in that it uses the Fredholm form of the 



integral equation (whose range of integration is from to rmax ), as done in Ref. 



and does not transform it into a Volterra type (whose range of integration is from to the 
variable radial distance r). It thus it avoids the need to evaluate the constants C which occur 
in the Volterra method, but instead it has to solve for a larger number of other constants. 
The latter arise by dividing the radial integration interval [0 Tmax] into partitions, and by 
expanding the solution in each partition into two independent functions which in turn are 
obtained by solving a local integral equation through the spectral expansion into Chebyshev 
Polynomials. There are twice as many such coefficients as there are partitions, and hence the 
matrix from which the coefficients are calculated is large, say 600 x 600. However, this matrix 
is sparse, and hence soluble economically. As is shown here, the semi-separable structure 
of the exchange nonlocality allows us to preserve the sparseness in the present case as well. 
If, however, the nonlocal potential is not of the semi-separable form our integral equation 
method still gives stable and accurate solutions, but then the biggest matrix involved is 
no longer sparse [^]. From the numerical point of view, the main difference of our S-IEM 
from other integral equation methods described above is that the radial mesh-points in the 
S-IEM are not equidistant, while those for the latter are. The numerical errors of the latter 
are of the finite difference type, and are given by a fixed power of the distance between 
mesh points, while in the S-IEM the errors become smaller than any power of the distance 
between mesh points. Or, more precisely, in the S-IEM the errors become smaller than any 
inverse power of the number of Chebyshev support points in each partition, a property which 
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expresses the spectral nature of the S-IEM. This property permits the S-IEM to have far 
fewer mesh points than the more conventional discretization formulations of differential or 
integral equations for a given envisaged accuracy, as will be demonstrated in Fig. ^ below. 

The first of our three comparison methods consists in replacing the exchange kernel by 
a small number of fully separable terms, and carrying out iterations only over the remainder. 
This is possible because, as is well known, the Green's function for a Schrodinger equation 
with both local and non local but fully separable potentials can be obtained without much 
difficulty by adding terms to the Green's function distorted only by the local potential. Our 
second method is a modified integral equation method, denoted as M-IEM. It uses a set of 
basis functions which are obtained by applying successively higher powers of the hamiltonian 
operator with local potentials on an initial scattering wave function. This method has been 
very successful in applications to nuclear physics problems. The third method the S- 
lEM, uses the Lippmann-Schwinger integral form of the Schrodinger equation. It differs 
from a previously introduced non-iterative solution of the Lippmann-Schwinger integral 
equation, denoted as NIEM |[T^, in that it uses non-equidistant mesh points, divides the 
radial interval into partitions of adjustable size, and uses a very accurate spectral integration 
technique involving Chebyshev polynomials. Because of its inherent stability the S-IEM is 
likely the method of choice for situations requiring solutions out to large distances. The 
generalization of this method so as to include the exchange potential [Q, [^] is described 
further below. 

In section 2 the basic equation to be solved will be described; in sections 3, 4 and 5 
the S-IEM, the SVD-improved iterative method, and the method of moments, respectively, 
will be reviewed; in section 6 the numerical comparison between the four methods will be 
described; section 7 contains the summary and conclusion; and Appendix 1 contains further 
details of the extension of the S-IEM to the presence of exchange. 



II. EQUATIONS AND NOTATIONS. 



The equation describing two electrons, one bound to a hydrogen-like nucleus of charge Z, 
and another incident with kinetic energy Ek on the ground state of the atom is 



ri 



+ — 

r2 



^(fi,r2) = E^(ri,r2), 



(5) 
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where \l/(ri,r2) is the overall wave function, E is the total energy, e is the charge of the 
electron, fi is the is the reduced mass of the incident electron, h is Planck's constant, fi and 
r2 denote the position vectors of the two electrons, ri and r2 are the respective magnitudes, 
and ri2 = \fi — is the distance between the two electrons. In order to transform the 
variables into atomic units, one multiplies Eq.(|^) by {2fi/h?) a^, where ao = {h'^/fie'^) is the 
Bohr unit of length, with the result 

^(fl,f2) = E^(fi,f2). (6) 

Here x = r/ao is a displacement vector in units of Bohr, and E = E/^ is the total energy 
in Rydberg units, with 3? = h'^/{2fial). 

In the Hartree-Fock approximation one expands the total wave function in terms of the 
bound states (pi, i = 1,2, ... of the atomic electron 

^(n, r2) = J2 <^i(^2) ± Mr2) Mn)] , (7) 

i 

where ipi are the wave functions of the scattered electron in channel z , to be determined from 
the solution of a set of coupled equations. The + or the — signs occur for the spin singlet or 
triplet cases, respectively. The subscript i represents the set of all quantum numbers which 
label the electron bound states. The corresponding principal quantum number is rii, and 
the corresponding bound state energy is = —{Z'^/nf)'^. The case of two or more bound 
electron states can also be derived. The result is a set of coupled equations with local and 
non-local pieces in the diagonal and off diagonal potentials. The latter are semi-separable 
of fully separable, hence the method described here for the one-channel case can also be 
applied. In the present study only the ground state will be assumed, i.e., i = 1 , and 
henceforth this subscript will be dropped, and further, Z = 1 . Under these assumptions 
the bound-state electron energy is e = — 3? and the incident electron has the asymptotic 
kinetic energy E^ = E — e . Assuming that this is a positive quantity, the corresponding 
wave number k in units of ao is given by 

e = Ek = (E- e^) /^ = E-l. (8) 

The equation for ip is obtained by truncating the sum in Eq. (|^) to one term, inserting 
it into Eq. (^, multiplying on the left by the functions (f) = 4'i{'T*2), and integrating over 



X2 



2Z 2Z ^ 2 

Xi X2 Xi2 
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d^r2. In the present numerical study only the case of orbital angular momenta = will be 
considered. The result is 



+ V{x,) - k^)^{n) ± + e)] <(j)\^^J> + < </)| — > = 0, 

3^12 



(9) 



< (P\^4> > and the symbol < | > denotes < A\B >= J A{r2)B{f2) d^r2. 



where V - 

If, furthermore, the radial wave functions Rl are introduced in the usual way 



ri ^-^ 



(10) 



where the the 's are Legendre Polynomials, and if Eq. (|) is multiplied by Pi^{co^Q\) and 
integrated over the solid angle d£lx one obtains the final equation for the radial function Rq 
for L = 



In the above, (assuming Z = 1), 



Roixi) = V{xi) Ro{xi) ± I ^{xi,X2) Roix2) dx2. 





V{x) = -2e^'^(l + 

S(xi,X2) = u{xi) U{X2) 

7 = -A;2 - 1, 

u{x) = 2a;e"^, 
, . 1 



X 

2 ■ 

7 + — 

Xl2 



vix) 



-u(x) = 2e-^. 



X 



(11) 

(12) 

(13) 

(14) 
(15) 
(16) 



The result for u arises from the well known expression for (pi 

01 = {Z/aof^"^ 2exp{-Zx2)Yoo{f2) 

with Z = 1. Utilizing the expansion of l/xi2 into Legendre Polynomials in the angle between 
the directions of xi and X2 , and remembering that only the term in Pq enters in the present 
case, one can recast the kernel in the semi-separable form 



55(xi, X2) = 2v{xi)u{x2) + ju{xi) u{x2) forx2 < Xi 
53(xi, X2) = 2u{xi)v{x2) + 7'u(xi) ^(^2) for X2 > xi. 



(17) 
(18) 



The above equations ([T7| ) and ([TSl) are the ones which will be used by the three methods 
of calculation, to be described in sections 3, 4, and 5 below. 



10 



It is interesting to note that Eq. (|TT]) can be replaced by an equivalent set of coupled 
equations[p^, pT 



+ e- V{x) 



dx"^ 



Ro{x) = ±Vi2{x)ip2{x) ±c^u{x) 
ip2{x) = V2i{x)Ro{x) 

u{x')RQ{x')dx' 



(19) 



with 



^12(2;) = ^21(3;) = V8exp(-x). 



(20) 



The reason that it is possible here to replace a nonlocality by an equivalent added channel 
is that the Green's function which corresponds to the operator d'^/dx'^ is given by the product 
/(x<)(7(x>) , with /(x) = X and g{x) = 1. This equivalence is due to the fact that the 
Coulomb interaction l/ri2, which appears in the first nonlocal term in Q", is closely related to 
the Laplacian d'^/dx'^. For angular momenta L other than zero it is sufficient to add the term 
L{L + l)/x^ into the square bracket of the second equation above. Whether the addition of 
extra channels is feasible for exchange interactions different from the Coulomb interaction 
remains to be investigated. 

The above equations (|19D are a set of inhomogeneous coupled equations which can be 
solved by conventional numerical means. However for more general semi-separable non- 
localities, which cannot be reduced to a set of equivalent coupled equations, the methods 
of solving Eq. ([TT|) presented in the next sections can be used. The question of whether 
an exchange nonlocality gives effects which are similar to a nonlocality due to coupling to 
inelastic channels has been examined by many authors. For certain nuclear scattering cases 



the two nonlocalities gave quite different results||28|. It is easy to understand the difference 
from Eqs. (|19D above, since in the exchange case the second channel contains no energy, 
and a inhomogeneous term is present, while both features are absent in the inelastic coupled 
channel case. 



11 



III. THE INTEGRAL EQUATION METHOD (S-IEM) 



In this section we describe how the previously developed spectral integral equation 
method for local potentials [Q can be extended so as to include the exchange terms. We 
drop the channel subscripts, since the equation being discussed contains only one channel, 
and for ease of notation we replace the function Ro{x) in equation (|TTD by if{x). In the 
integral equation method S-IEM the differential equation (|Tl|) is transformed into the 
equivalent integral equation 

/»00 /"OO 

(p{x) = sm{kx) + / Q{x, x')V{x')(p{x')dx' + / ^{x, x")(p{x")dx" . (21) 
Jo Jo 

where Q is the undistorted Green's function corresponding to the momentum k, 

Q{x,x') = ——cos{kx) sin(A;x') for x' < x 

Q{x,x') = ——sm{kx) cos{kx') for x' > x, (22) 

rC 

and the kernel results from the presence of the exchange terms, and is equal to the 
convolution of the Green's function with the nonlocal potential Q", defined in Eqs. ([T7|) and 

POO 

J-'{x,x")= / Q{x, x')'^{x' , x") dx' . (23) 
Jo 

For a general nonlocal potential the kernel JF is not semi-separable. In this latter case 
the solution of the integral equation can still be performed and gives rise to matrices which, 
although not sparse, have a structure such that they can still be evaluated economically 
26| . However, if the nonlocality 53 is semi-separable, as is the case when it results from 



exchange terms, then it can be shown |]5[, that the kernel , Eq. (p3D also is semi- 
separable and is of rank 2. Even though the resulting expression for JF is not as simple as 



the rank 1 expression (p2D for Q , the conventional lEM method for local potentials can be 
extended to this will be shown below. The advantage of this technique is that the 

"big" matrix, which occurs in the process of piecing together the local solutions obtained for 
each partition, is a sparse band limited matrix, and hence the complexity of the calculation 
remains proportional to the number of partitions m, rather than being of power m^, as would 
be the case with general non-sparse matrices. However, the complexity of the calculation 
also contains a factor which increases like the cube of the number of bands, which, in the 

12 



case of the presence of nonlocal exchange potential doubles compared to the local case, and 
hence the complexity of the calculation increases by an overall factor of eight. 

The semi-separable form of the kernel T is obtained by inserting into Eq. ( p3D the 
expression ( P^j ) for the Green's function, and using for ^{x\ x") the expressions given by 
Eqs. and If one also combines the integral over QV in Eq. (^) with the kernel 

JF into a single kernel K 

POO 

ip{x) = sm{kx) + / K{x, x")^{x")dx\ (24) 
Jo 

one obtains for K the result 

K{x, x") = h{x) gi{x") + f2{x) g2{x") for x" < x (25) 
K{x, x") = pi{x) qi {x") + p2 (x) q2 {x") for x" > x. (26) 

The subscript 1 and 2 stands for the first and second semi-separable terms in the rank-2 
expression for K , respectively. The functions /, p and q are given by 

/i(x) = cos(A;x) (27) 

gi{x) = 2Xsu{x) v{x) — 2Isv{x) u{x) — — sin(/cx) V{x) (28) 

f^[x) = cos{kx) [2Isy{x) + 7Xsn(x)] + sin(A;x) [2Icv{x) + jTcu{x)] (29) 
g2ix)=u{x) (30) 

Pi{x) = sin(fca;) (31) 
qi{x) = 2Xcv{x) u{x) — 2Xcu{x) v{x) — — cos{kx) V{x) (32) 

P2{x) = cos{kx)2su{x) + sm{kx)2cu{x) (33) 
q2{x) = 2v{x) + 7m(x), (34) 

where the functions T are defined by 

1 r 

Isuix) = —-r I sm.ikr)u{r)dr 
k Jo 



1 r 

^svix) = —-r / sm(kr)v(r)dr 
k Jo 



k 



oo 



1cu{x) = —— I cos{kr) u{r) dr 



X 

oo 



1cv{x) = —J I cos{kr) v{r) dr (35) 
k Jx 

and the functions u and v are defined in Eqs. (p^Sj) and ([T6|), respectively. 
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A. The discretization. 



The integral equation to be solved is 



{x) + ^^ r gi{x')<P{x')dx' + I g2{x')(P{x')dx' 



^ Jo ^ Jo 

+^4^ r (liix')H^')d^' + r q2{x')(l){x')dx' = sin{kx). 

It can be written in concise form as 

{I + K)(l){x) = sin{kx). (36) 

The radial distance x is contained in the range < x < rmax ? where the upper limit r^ax is 
chosen sufficiently large so that beyond rmax the integrands can be neglected. We start by 
dividing the interval [0, rmax] into m partitions [bo, bi], 62], • • • , [k-i, k], ... [6^-1, &m], 
where the points bi are not necessarily equispaced, similarly to what was done in Next 
we show that in each interval i the global solution (p can be found as a linear combination of 
four local solutions of equation (|36|) restricted to each of the subintervals of partition. Let 
Ki denote the operator K restricted to act only in the subinterval bi]. For example, Ki 
operating on the function 77 is given by 

{K,ri){x) = ^^ r g,{x')ri{x')dx' + f\i{x')7]{x')dx' 



+ ^ r 92{x')Tl{x')dx' + r q2{x'Hx')dx', 



bi-i <x<bi. 

Then, in terms of Ki, equation ( pG]) can be rewritten as 

(/ + Ki)(P{x) = + 5«p2(x) + + /^«/2(x), (37) 

bi-i < X <bi, 

where use has been made of the fact that pi{x) = sin{kx). This result can be obtained 
(see Ref. [||) by decomposing the integrals in equation (^) into three domains: [0, 

and rmax]-The second domain gives rise to the operator K^. Accordingly the 
constants are given by 

I A" ''max 

A« = 1 - 7 / qiix')(f){x')dx', (38) 
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k 



q2{x')(f>{x')dx' , 



gi{x')(f){x')dx' , 



(39) 
(40) 



g2ix')(Pix')dx', (41) 

which are later found from the solution of matrix equation (^) We next define four functions 
yi,Zi, fii, in each sub interval i by 

{I + K,)y,{x) =pi{x), (42) 
{I + K,)z,{x) = fi{x), (43) 
{I + Ki)fii{x) =p2{x), (44) 
iI + K,)Ux) = f2{x). (45) 

In view of the fact that the operator Ki is linear, the solution 0(x) of equation ( ]37| ) in each 
subinterval i is given by 



< X <bi. 



(46) 



This result allows one to relate the constants A, B, C, D in subinterval i with those in other 



subintervals j, by inserting ( ^61 ) into equations (p8D-(pD. The resulting equations, described 
in Appendix 1, can be transformed into a block tridiagonal-system, similarly to what was 
done in our previous work, 



I 


U12 














" Ai" 




" 0' 


U21 


I 


U23 











A2 










U32 


I 


U34 









A3 















I 


^ m—l,ir>. 



















Um,m— 1 


I 




Am 




E 



(47) 



where the quantities A, 0, and are 1x4 column vectors 

A, = [A«,5(*),C(*\d(*)]^, 
0= [0,0,0,0]^, 
E= [1,0,0,0]^ 
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and each block is a 4 x 4 matrix, 



U, 



(a2l/)i+i («2/u)i+i - 1 {a2z)i+i (a20i+i 






(48) 



u 



i+l,i 



(49) 







(A/i). {PlZ)^ - 1 (AO. 

for i = l,...,m— 1, and I is a 4 x 4 identity matrix. The entries of the matrices are integrals 

in each partition i of products of the known functions qi,q2, gi, g2, and the numerically 
computed functions y, fi, z, and ^ , where by definition, 

-6,; 



p{x) q{x) dx. 



It is noteworthy that the structure of the matrix in Eq.(^7D is very similar to the structure 
encountered for a set of four coupled channels (see Eq. (34) in Ref. [0]. (This reference 
contains further details of the discretization technique). The system of equations (^71) is 
solved by Gaussian elimination specialized for band limited matrices, see e.g. Its 
complexity is 4mp(p + 1) — |j9^, where p is the number of non-zero subdiagonals, (band- 



width), and m is the number of partitions. In Eq. (|47|) p = 7. Since the number of grid 
points per partition, 16, is larger than p = 7, it is clear that the overall cost will be dominated 
by the cost of solving Eqs. (^2]) in all partitions, which is of order 16^m 



IV. THE SVD-IMPROVED ITERATIVE METHOD. 



The first of our three comparison methods consists in replacing the exchange kernel by a 
small number of fully separable terms, and carrying out iterations only over the remainder, 



as will be described in this section, and as is given with more detail in Ref. |T^. As is 
well known, the Green's function for a Schrodinger equation with both local and non local 
but fully separable potentials can be obtained without much difficulty by adding terms to 
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the Green's function distorted only by the local potential. By contrast, the older iterative 
method of taking the nonlocal kernel into account perturbatively consists in writing Eq. 
(pAl) in the form 



dx{ 



2 -V{xi) + k'' 



Rq{Xi)=± / ^{Xi,X2) Rq{x2) dX2 



(50) 



and then transforming it into the iterative integral equation 



/■oo / /"OO \ 

/2o"^^^(a;i) = /(xi) + j Qv{xi,x) \ ± J '^{x , X2) R^q\x2) dx2j dx . 



(51) 



In the above, /(x) is the "regular" solution of 



dx{ 



2 -V{xi) + k'' 



f{x)=0, 



(52) 



and ^y(xi,x') is the Green's function which corresponds to the left hand side of Eq. (pOD. 
It is distorted by the local potential V, and can be expressed in terms of semi-separable 
expressions involving two independent solutions /(x) and g{x) of Eq. (|52D, 

1 



Gv{x,x') 



f{x)g{x') for X < x' 



Qv{x,x') = —— g{x)f{x') for x > x', 



(53) 
(54) 



as is well known. The functions / and g are normalized such that their Wronskian is equal 
to k. The iteration is started by using the solution in the absence of the exchange terms for 
the first {n = 0) guess R*^^ = /(X2). 

The rate of convergence of the iterations depends on the norm of 



POO 

J-'v{x,x")= / Gv{x, x')^{x' , x") dx' . 
Jo 



This norm in turn depends on the norm of 3", and on the norm of Qv The latter becomes 
large at small incident energies fc^ , in view of the presence of the factor 1/ A; in Eq. (|53|) , and 
hence the iteration will diverge for a sufficiently small value of k. The rate of convergence 
also depends on the ± sign in front of the exchange integrals, as was found in the numerical 
examples described below. This effect does not occur in the other methods described in this 
paper because the latter do not make use of the iteration on T^v- 

In what follows in this section, we describe a method, to be denoted as SVD, which 
reduces the norm of the nonlocal kernel by decomposing it into a sum of a fully separable 
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kernel of low rank plus a remainder. The separable terms are placed in the left hand side of 
Eq. (0), the Green's function in the presence of both the local distorting potential V and 
the separable nonlocal pieces of the kernel is obtained, and hence iterations of the form of 
Eq. (^Tj) can be carried out, where $5 is now the residual kernel. This way the "/c" -divergence 
can be shifted to smaller values of k, but it cannot be avoided. 

A. The separable content of the nonlocal kernel. 

The singular value decomposition method (SVD) is used to decompose the kernel 
Q{xi,X2) into a number of fully separable terms plus a remainder. The method is as fol- 
lows. First a numerical integration algorithm is chosen which divides the range of integration 
[0, -Rmax] into a set of discrete points. Correspondingly the kernel X2) is transformed 
into a N X N matrix K{i,j), with i,j = 1, 2, ..A^. We next perform a singular value decom- 
position on K. The SVD method is based on a theorem of linear algebra according to which 
any M x N matrix K can be written as the product of an M x M orthogonal matrix U, 
an M X N diagonal matrix S with positive or zero elements, and the transpose of an 
X orthogonal matrix V. (A matrix U is orthogonal if UU^ = U'^U = I, which means 
that its columns are normalized and orthogonal to each other, and so are the rows.) For our 
purpose it is sufficient to consider the case A^ = M. In this case we can rigorously write 

TV 

K = U^V^ = J2 (^sUs^rJ (55) 

where the columns of U and V are the column vectors Ug, and v^, respectively, and S is a 
diagonal matrix of the non-negative quantities as, s = 1,2, ...N, ordered by decreasing size 
(the largest ones first). The latter are the "singular values". As a result of the above, a fully 
separable piece of rank n can be separated out of the matrix K, leaving a residual matrix 

K = + K^. (56) 
by carrying the sum in Eq.(|55|) to a upper limit n which includes only the largest values as- 

n n n 

K^ihi) = ^Ujs(Ts Vsi , or K-'^ = ^ o-^Ugvf = ^Ug) as (v^. (57) 

s=l s=l s=l 
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The last entry into the above equation uses the Dirac notation for a vector and its transpose. 
The remainder is given by 

N 

K'' = K-K^ = ^sUs vf . (58) 

s=n+l 



B. Greens function for a separable potential. 



In order to obtain the Green's function Qv+k={x, x'), which is distorted by both the local 
potential V and the fully separable Kernel , we rewrite Eq. (p^) symbolically in the form 



(59) 



where the integration over the variables is implicitly assumed. For simplicity, let us assume 
that only two terms in are responsible for the divergence of the iterative Green's function 
approach, Eq. (^Tp. In order to obtain the overlap integrals = 1,2 we multiply 

Eq. ( |5DD on the left with y/ai{vi and integrate over all x's, with the result that y/oKviip) = 
■\f^i{vif)+^/ai{viQv{K^ . Rearranging terms one obtains the following matrix equation 
for JaU^Viii) 





= \fo 


' ivif) ' 

















(60) 



where 



M 



and 



1 + 6^11 


Ql2 1 


















6^21 1 


+ Q22 ) 








Gij = \f^i 


{ViGv u 




= 1,2. 





Solving Eq. (|60| ) for [{viip), {v2ip)] and inserting the result into Eq. (|59D , one obtains 
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from which the result for Qv+k= emerges: 



1 - 



Ml) U2) 



(^1 

{V2 



(62) 



The numerical result for the triplet phase shift, shown in Table 1 of section 6, used five 
sets of singular value functions and and required five iterations of Eq. |51|. Without 
the use of the SVD expansion, the iterations did not converge for k < 0.3. For values of 
k < 0.1 the iterations using the SVD expansion did not converge for either the singlet or 
triplet cases. 

V. THE MODIFIED INTEGRAL EQUATION METHOD (M-IEM) 

In this section we describe the method proposed by Kim and Udagawa[^, which we call 
the modified integral equation method (M-IEM). The method is well documented in the 
literature, and hence only a brief description is given here. It starts from the following 
equation obtained by rewriting Eq. fpUD ; 



ip{x) = ±A(x) 



A(x) = / ^{x, x') (p{x') dx' . 
Jo 

We then transform the equation into the integral form as 

ifix) = ip^^\x) ± / g'{x,x")\{x")dx", 



where (f^^\x) and Q'{x,x') satisfy 



d 



dx\ 



dx{ 



2 -V{x) + k'' 



¥^(°)(x) = 
g\x,x") = 6{x-x") 



(63) 
(64) 

(65) 

(66) 
(67) 



Further, we modify Eq.(^) by multiplying both sides by 53(x,x') and carrying out the 
integration over x' . The result is 



00 fOO 



\{x) = X^^\x)± / '^{x,x')G'{x\x")\{x")dx"dx' 
Jo Jo 



AW( 



X] = I ^{x,x')ip^°\x')dx' 





(68) 
(69) 
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The equation we solve is (|6^) . Since both ip^^\x) and Q'(x,x") are defined in terms of the 
local potential V{x), they can be calculated without any problem. This means that once 
the solution A(x) of Eg . (|68D is obtained, then ip{x) can be calculated from Eq. (|65D . 

In solving Eq. ([68|) , use is made of the Lanczos method ||31|. It is worth noting that the 
application of the Lanczos method for solving Eq.(|68D is possible, since A(x) is a bounded 
function, as can be seen from the fact that it is essentially given in terms of the bounded 
nonlocal potential function Q{x,x'). This makes it possible to expand X{x) in terms of an 
orthonormal set of functions, as is done in Eq.(|74[) below. This is not the case for <f{x) in 
Eq. (|65D , since (p{x) is not bounded. 

We first expand X{x) in terms of the orthonormal set of functions Di{x) with i = 
0, 1, 2, ,Ni which are generated as follows: 



"0 

1 j;^ j;^'^{x,x')g'{x',x")D^^i{x")dx"dx'- 



di 



(70) 
(71) 



with 



/o" /o" /o" D,{x)%x, x')g'{x', x")D,{x") dx"dx'dx, j <i + l 



(72) 



j >i + l 

The normalization constant d^ in Eqs. (ffOD and (|7l| ) is determined from the condition 

D,{x)Di{x)dx = 1 (73) 



Di{x) being the conjugate function to Di{x). The coefficients aji given by Eq.(^) are those 
determined from the usual Schmidt orthonormalization procedure. Now we write X{x) as 

Ni 

X{x) = J2c,D,{x), (74) 

j=0 

where Cj are the expansion coefficients. 

Inserting Eq. ([7^) into Eq. (|68|) , one can easily derives a set of inhomogeneous linear equa- 
tions for the expansion coefficients Cj, i.e.. 



^^(5jj — aij)Cj — doSoi. 



(75) 
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0.0 0.2 0.4 0.6 0.8 1.0 

FIG. 1: Dependence of the singlet (+) and triplet (-) L = phase shifts on the incident wave 
number, displayed in the form k/tan{6). Where these curves pass through zero, tan((5) goes 
through infinity, as is shown in the next two figures. Neither the M-IEM or the S-IEM had any 
difficulty evaluating these quantities either for the small values of k or in the vicinity of the zeros. 



The values of Cj are then determined by solving Eg. ([75|) . Note that Eq.(^) can be solved 
rather easily, because aji = for j > i + 1 (see Eq.(|72D). In addition, the value of Ni can be 
chosen as a small number. This helps greatly in making the actual numerical calculations 
very fast. 



VI. NUMERICAL RESULTS 

The bench-test calculation performed by the three methods described above consists in 
obtaining the L = phase shift for the scattering of an electron from the ground-state of an 
Hydrogen atom, in the presence of exchange terms, both for the singlet and the triplet states, 
S^^^ and S^~^ , respectively. The methods are the S-IEM, the SVD, and the M-IEM. The 
older integral equation method is denoted as NIEM (the "N" stands for non-iterative), and 
a representative result is taken from the paper by Sams and Kouri [§, since these authors 
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0.24 0.26 0.27 0.28 0.30 



FIG. 2: Momentum dependence of the singlet phase shift in the vicinity of 7r/2 (mod. vr), as 
calculated by the S-IEM. Please note the scale of the y-axis. 



describe their accuracy for exactly the same test case as ours. The SVD, M-IEM and the 
NIEM methods use equi-spaced mesh points, since their auxiliary functions are the solutions 
of local differential equations using finite difference methods, while the S-IEM, as mentioned 
above, uses non-equispaced mesh points, which are the zeros of a Chebyshev polynomial of 
a certain order (16 in this case), in each of the partitions into which the radial interval is 
decomposed. 

The k dependence of the phase shifts is shown in Fig. 1, by plotting the ratio k/ tan(5). 
When either of the two curves crosses the line, the corresponding value of tan(5) becomes 
infinite, as is shown in Figs. (0) and (|^), and the respective phases shift have the value tt/2, 
modulus vr. Both the M-IEM and the S-IEM methods had no difficulty in reproducing the 
singularity in tan(5), and both were able to reach arbitrarily small values of the momentum 
k. By contrast, the SVD method could not obtain results for k < 1.0(ao)~^. 

Since well documented accuracy studies exist for the NIEM p we examined the rate of 
convergence of the phase shift as a function of number the mesh points for a case which 
is treated in Ref. [§. The case chosen is the singlet phase shift with exchange, S^'^\ The 
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FIG. 3: Same as Fig. 2 for the triplet phase shift. 



value of the wave number is k = 0.2 (ao)~^ , and the maximum radial distance is 20 oq. The 
number of significant figures obtained for each of the non-spectral methods, all using a mesh 
size of 0.005 oq, is shown in the three last rows of Table 1. The result for the spectral method 
S-IEM, also shown. 

Table 1: Accuracy of S^'^^ for various algorithms. 



Method 




# of Pts. 


S-IEM 


1.8701579 


80 


M-IEM'^) 


1.870156 


4000 


NIEM ^) 


1.87015 


4000 


SVD 


1.8701 


4000 



"■^ Five basis states D are used in this calculation 
^^Non- iterative method of Ref. 



The convergence of the four methods with the number of mesh-points is illustrated in 
Fig. ^. The number of significant figures for a given number of mesh points is determined 
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from the stability of the result obtained after rounding, when compared to the result with 
the next higher number of points. It is clear from the figure that the S-IEM method reaches 
higher accuracy with a smaller number of points than the other methods shown. With 160 
mesh points (the corresponding number of partitions is 10) the value obtained for 5'^^^ = 
1.87015788462442 rad. is the same , to within the quoted number of 15 significant figures, 
as the result for 224 mesh points. This is close to machine accuracy, and shows that the 
accumulation of round-off errors is small in the S-IEM method, confirming previous studies. 
(See Fig. 1 in the 1997 paper quoted in Ref. [Q]). The discrepancy in the seventh significant 
figure between the S-IEM and the M-IEM could be due to the fact that only five basis 
functions were used for the latter. This point has not been investigated further. 

The scattering length a and effective range for this one state electron-hydrogen scat- 
tering calculation have also been examined. The procedure is similar to the one used for a 
previous atom-atom scattering bench-mark calculation 0. It is based on the low momentum 
expansion of the scattering phase shift 

kcot6o = + reP + 0{k^). (76) 
a 

The left hand side of the above expression is calculated for two very small values of the 
wave number k , differing by a factor of two, and the values of a and are solved for. The 
procedure is repeated for decreasing values of k and increasing values of the maximum radial 
distance rmax and of the number of mesh-points until stability in the results is found to a 
given number of significant figures. For the values of a and rg listed in the tables below for 
the S-IEM method, values of A; ~ 10~^ , rmax ~ 50 and approximately 1000 mesh points 
were found to be adequate. However, contrary to what was done in Ref. the value of 
'"max was not extrapolated to oo via a perturbative method. For the M-IEM case, the values 
of a and Tg in Table IV are extracted from Eq. |76| by calculating k cot 6q for the two values 
of /c = 0.00001 and 0.01, and for T = 20. Excellent agreement with the S-IEM values is 
obtained. 



Table II: Scattering lengths a. 



Method 


Singlet 


No exchange 


Triplet 


S-IEM 


8.100312397 


-9.44716668854 


2.349396156 


M-IEM 


8.1003 


-9.44716 


2.3494 
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FIG. 4: Comparison of the accuracy of the S-IEM, M-IEM and NIEM integral equation methods 
for the calculation of the singlet phase shift, as a function of the number of meshpoint used. The 
incident momentum is A; = 0.2(ao)~^, and the value of the radial cut-off point is rmax = 20 ao- The 
NIEM results are taken from Ref. The accuracy for a given number of mesh points of each 
method is determined by the number of significant figures which are stable (after rounding), as 
compared with the result for the next higher number of meshpoints. 



Table III: Effective Range r^. 



Method 


Singlet 


No exch. 


Triplet 


S-IEM 


1.51201 


0.766797 


0.6105 


M-IEM 


1.51 


0.767 


0.612 



VII. SUMMARY AND CONCLUSIONS. 

In this paper four methods were compared to solve the one-dimensional Schrodinger 
equation in the presence of the exchange nonlocality for the case of electron scattering from 
a hydrogen atom, with only the lowest energy state of the bound electron being included. 
The oldest method in the literature proceeds by first solving the equation in the presence 
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of only the local potential, and then including the nonlocal part through Green's function 
iteration. The iterations converge only for a limited range of parameters, and one of the 
methods described here improves upon the convergence by separating out of the nonlocal 
kernel a fully separable part by means of the Singular Value Decomposition method (SVD). 
By this means the region of convergence could be extended to a larger domain, but for our 
example convergence still fails at small values of the momentum, k < 0.1(ao)~^ and the 
maximum accuracy achieved was five significant figures. Another method was developed 
in the literature, in which the differential Schrodinger equation is first transformed into 
an Lippman-Schwinger integral equation, and is then solved non-iteratively (NIEM). In 
our example, taken from the literature 0, this method achieved six significant figures of 
accuracy, but appears not to work for small values of the incident momentum. Improved 
accuracy and the viability for all values of k was achieved in the present study by extending 
a previously developed spectral solution of a Lippman-Schwinger integral equation with 
local potentials [|l|, [@, to the case with an exchange-type nonlocality. This extension was 
possible because the exchange nonlocality is of a semi-separable character. The resulting 
method (S-IEM) gives substantially higher accuracy (15 significant figures) than the NIEM, 
and converges much faster with the number of mesh-point in the integration interval than 
the NIEM, as is illustrated in Fig. 4. A fourth method, (M-IEM) developed previously 
for nonlocalities occurring in nuclear physics ||^ achieves seven figures of accuracy, and has 
no difficulty in coping with small values of k. The rate of convergence of this method was 
comparable to that of the NIEM. The reason is due to the fact that the auxiliary functions 
needed for both methods, as well as the integration algorithms, are based on a finite difference 
algorithm, whose error usually decreases inversely with the number of mesh points according 
to a well defined power. From inspection of Fig. 4, this power has the relatively low value of 
2.5. Both the M-IEM and the SVD methods have the advantage that they can be used for 
non-localities which are more general than the semi-separable exchange ones. The S-IEM 
also can be applied to these cases, but, at its present stage of development, the large matrix 



in Eq. ^ is then no longer sparse p6 



In summary, four methods of calculating the scattering phase shift in electron atom 
collision were compared for a numerical test case, and the advantages and disadvantages of 
each were discussed. 
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Appendix 1 

The relation between all coefficients in fH6D can be written in the matrix form, 





M12 


Mi3 


Mm 




A 




T 


M21 


M22 


M23 


M24 




B 







M31 


M32 


M33 


M34 




C 







M41 


M42 


M43 


M44 




D 








eg : 15 



(77) 



where 



A= B=[B,,...,B, 



T 



T 

mj ) 



C = [Ci,...,CX D = [Di,...,D, 

i = [i,..,if o = [o,...,o]^, 



and each of the block matrices Mij are either upper triangular matrices U or lower triangular 
matrices L with either I's or O's in the main diagonal, respectively (the subscripts are 
accordingly 1 or ). These matrices Mj are thus of the form 



1 72 72 ■ ■ ■ 7n 



1 73 



7n 



■■■ 1 7„ 
••• 1 



f/n 



72 72 
73 



7m 
7m 



■■■ 7^ 
■■■ 



or 



1 ■■■ 

<5i 1 ■•• 

6i 62 1 

61 62 ■ ■■ 



5. 






■■ 

1 

n-1 1 



■■■ 

61 ••• 

Si 62 

61 62 ■■ ■ 






... 



6m-l 



in which the entries 7 or 5 are given in the table below 

Table 1: Entries 7 and 6 in the block matrices M=U or M = L. 
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A/ii = f/i with = {oav^^ 


ALio = f/n with = foi u)i 


Afi*? = f/n with 'y^ = (o-\z)^ 


/Vfi/i = f/n with = foi^^li 


Afoi = r/n with = (ooii^^ 


/VLoo = r/i with o/,- = iooii?\i 


Mori = (In with 'y^ = (ooz)^ 


jA/fo/i = IJn with = iO'^f\i 


M'i^ = Ln with 5.- = foiw),- 


= Ln with 5i = folu),- 


M33 = Li with (5i = {giz)i 


M34 = Lo with (5i = {giCji 


M41 = Lo with 5i = {g2y)i 


M42 = Lo with 5i = (5'2/i)j 


M43 = Lo with 5i = (5(2 2:) i 


M44 = Li with 5i = (fifaOi 



Using elementary row operations on equation (^) and then changing the order of the vari- 



ables, the coefficient matrix of equation ([77|) can be transformed into the block tridiagonal 



system (|T7|) to (|49|), given in the text. Further details can be found in Ref. [|I| 
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